import numpy as np
import pandas as pd
import geopandas as gpd
import math
import folium
from folium import Choropleth, Circle, Marker, GeoJson
from folium.plugins import HeatMap, MarkerCluster
from geopy.geocoders import Nominatim
import warnings
warnings.filterwarnings('ignore')
import os
print(os.getcwd())g:\내 드라이브\202301\Blog\Data_Mining\Geospatial_Analysis이 과정에서는 지리 공간 데이터 또는 지리적 위치가 있는 데이터를 다루고 시각화하는 다양한 방법에 대해 알아봅니다.
그 과정에서 다음과 같은 몇 가지 실제 문제에 대한 솔루션을 제공하게 됩니다.
또한, 보스턴의 범죄를 시각화하고, 가나의 의료 시설을 조사하고, 유럽의 최고 대학을 탐색하고, 미국의 독성 화학물질 방출을 추적할 수 있습니다.
이 첫번째 튜토리얼에서는 이 강좌를 완료하는데 필요한 전제 조건을 빠르게 다룹니다. 더 깊이 잇는 복습을 원하신다면, Pandas 강좌를 추천합니다.
이제, 첫번째 지리공간 데이터 집합을 시각화해보겠습니다.
첫번째 단계는 지리공간 데이터를 불러오는 것입니다. 이를 위해 GeoPandas 라이브러리를 사용하겠습니다.
지리공간 파일 형식에는 shapefile, GeoJSON, KML, GPKG 등의 다양한 형식이 있습니다. 이 강좌에서는 그 차이점에 대해서는 다루지 않겠지만, 중요한 것들은 다음과 같습니다. - shapefile이 가장 흔한 파일 형식입니다. - 지리공간 파일의 모든 형식은 gpd.read_file() 함수를 이용하여 빠르게 불러올 수 있습니다.
다음 셀은 뉴욕 주 환경보전국에서 관리하는 숲, 야생지대 및 기타 토지에 대한 정보가 포함된 shapefile을 불러옵니다.
# Read in the data
full_data = gpd.read_file("./geospatial-learn-course-data/DEC_lands/DEC_lands/DEC_lands.shp")
# View the first five rows of the data
full_data.head()| OBJECTID | CATEGORY | UNIT | FACILITY | CLASS | UMP | DESCRIPTIO | REGION | COUNTY | URL | SOURCE | UPDATE_ | OFFICE | ACRES | LANDS_UID | GREENCERT | SHAPE_AREA | SHAPE_LEN | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | FOR PRES DET PAR | CFP | HANCOCK FP DETACHED PARCEL | WILD FOREST | None | DELAWARE COUNTY DETACHED PARCEL | 4 | DELAWARE | http://www.dec.ny.gov/ | DELAWARE RPP | 5/12 | STAMFORD | 738.620192 | 103 | N | 2.990365e+06 | 7927.662385 | POLYGON ((486093.245 4635308.586, 486787.235 4... |
| 1 | 2 | FOR PRES DET PAR | CFP | HANCOCK FP DETACHED PARCEL | WILD FOREST | None | DELAWARE COUNTY DETACHED PARCEL | 4 | DELAWARE | http://www.dec.ny.gov/ | DELAWARE RPP | 5/12 | STAMFORD | 282.553140 | 1218 | N | 1.143940e+06 | 4776.375600 | POLYGON ((491931.514 4637416.256, 491305.424 4... |
| 2 | 3 | FOR PRES DET PAR | CFP | HANCOCK FP DETACHED PARCEL | WILD FOREST | None | DELAWARE COUNTY DETACHED PARCEL | 4 | DELAWARE | http://www.dec.ny.gov/ | DELAWARE RPP | 5/12 | STAMFORD | 234.291262 | 1780 | N | 9.485476e+05 | 5783.070364 | POLYGON ((486000.287 4635834.453, 485007.550 4... |
| 3 | 4 | FOR PRES DET PAR | CFP | GREENE COUNTY FP DETACHED PARCEL | WILD FOREST | None | None | 4 | GREENE | http://www.dec.ny.gov/ | GREENE RPP | 5/12 | STAMFORD | 450.106464 | 2060 | N | 1.822293e+06 | 7021.644833 | POLYGON ((541716.775 4675243.268, 541217.579 4... |
| 4 | 6 | FOREST PRESERVE | AFP | SARANAC LAKES WILD FOREST | WILD FOREST | SARANAC LAKES | None | 5 | ESSEX | http://www.dec.ny.gov/lands/22593.html | DECRP, ESSEX RPP | 12/96 | RAY BROOK | 69.702387 | 1517 | N | 2.821959e+05 | 2663.909932 | POLYGON ((583896.043 4909643.187, 583891.200 4... |
CLASS 열에서 볼 수 있듯이, 처음 5개 행은 각각 다른 숲에 해당합니다.
이 튜토리얼의 나머지 부분에서는 이 데이터를 사용하여 주말 캠핑 여행을 계획하는 시나리오를 고려해보겠습니다.
온라인에서 크라우드 소싱된 리뷰에 의존하는 대신 자신만의 지도를 만들어봅시다.
이렇게 하면 특정 관심사에 맞게 여행을 조정할 수 있습니다.
데이터의 처음 5개의 행을 보기 위해 head() 함수를 사용했습니다. 이 함수는 Pandas DataFrame에서도 사용된다는 것을 알고 있을 것입니다.
사실, 모든 명령어는 Pandas DataFrame에서도 함께 작동합니다.
이는, 데이터가 (Pandas) DataFrame의 모든 기능을 갖춘 (GeoPandas) GeoDataFrame 객체로 로드되었기 때문입니다.
type(full_data)geopandas.geodataframe.GeoDataFrame예를 들어, 모든 열을 사용하지 않으려는 경우 열의 하위 집합을 선택할 수 있습니다. (데이터를 선택하는 다른 방법을 복습하려면, Pandas 강좌의 튜토리얼을 확인하세요.)
data = full_data.loc[:, ["CLASS", "COUNTY", "geometry"]].copy()value_counts() 함수를 사용하여 다양한 토지 유형 목록과 함께 데이터 집합에 나타나는 횟수를 확인합니다. (이 함수 (관련 함수)를 복습하려면, Pandas 강좌의 튜토리얼을 확인하세요.)
# 각 유형별로 토지의 개수는 몇개인가요?
data['CLASS'].value_counts()WILD FOREST 965
INTENSIVE USE 108
PRIMITIVE 60
WILDERNESS 52
ADMINISTRATIVE 17
UNCLASSIFIED 7
HISTORIC 5
PRIMITIVE BICYCLE CORRIDOR 4
CANOE AREA 1
Name: CLASS, dtype: int64loc(iloc)를 사용하여 데이터의 하위 집합을 선택할 수도 있습니다. (이 함수를 복습하려면, Pandas 강좌의 튜토리얼을 확인하세요.)
# "WILD FOREST" 또는 "WILDERNESS"에 해당하는 토지를 선택합니다.
wild_lands = data.loc[data['CLASS'].isin(['WILD FOREST', 'WILDERNESS'])].copy()
wild_lands.head()| CLASS | COUNTY | geometry | |
|---|---|---|---|
| 0 | WILD FOREST | DELAWARE | POLYGON ((486093.245 4635308.586, 486787.235 4... |
| 1 | WILD FOREST | DELAWARE | POLYGON ((491931.514 4637416.256, 491305.424 4... |
| 2 | WILD FOREST | DELAWARE | POLYGON ((486000.287 4635834.453, 485007.550 4... |
| 3 | WILD FOREST | GREENE | POLYGON ((541716.775 4675243.268, 541217.579 4... |
| 4 | WILD FOREST | ESSEX | POLYGON ((583896.043 4909643.187, 583891.200 4... |
plot() 함수를 사용하여 데이터를 빠르게 시각화할 수 있습니다.
wild_lands.plot()<Axes: >
모든 GeoDataFrame에는 특별한 geometry 열이 포함되어 있습니다. 이 열에는 plot() 메서드를 호출할 때 표시되는 모든 기하학적 개체가 포함됩니다.
wild_lands['geometry'].head()0 POLYGON ((486093.245 4635308.586, 486787.235 4...
1 POLYGON ((491931.514 4637416.256, 491305.424 4...
2 POLYGON ((486000.287 4635834.453, 485007.550 4...
3 POLYGON ((541716.775 4675243.268, 541217.579 4...
4 POLYGON ((583896.043 4909643.187, 583891.200 4...
Name: geometry, dtype: geometry이 열에는 다양한 데이터 유형이 포함될 수 있지만, 일반적으로 Point, LineString, or Polygon 입니다.
이 데이터의 geometry 열에는 2,983개의 서로 다른 Polygon 객체가 포함되어 있고, 각 객체는 plot에서 서로 다른 모양에 해당합니다.
아래 셀을 통해 우리는 캠프장 위치 (Point), 도보 경로 (LineString), 지역 경계 (Polygon)를 포함한 3개의 GeoDataFrame을 더 생성합니다.
# 뉴욕 주의 캠프장 위치 (Point)
POI_data = gpd.read_file('./geospatial-learn-course-data/DEC_pointsinterest/DEC_pointsinterest/Decptsofinterest.shp')
campsites = POI_data.loc[POI_data['ASSET'] == 'PRIMITIVE CAMPSITE'].copy()
# 뉴욕 주의 도보 경로 (LineString)
roads_trails = gpd.read_file('./geospatial-learn-course-data/DEC_roadstrails/DEC_roadstrails/Decroadstrails.shp')
trails = roads_trails.loc[roads_trails['ASSET'] == 'FOOT TRAIL'].copy()
# 뉴욕 주의 지역 경계 (Polygon)
counties = gpd.read_file('./geospatial-learn-course-data/NY_county_boundaries/NY_county_boundaries/NY_county_boundaries.shp')다음으로, 4개의 GeoDataFrame을 모두 사용하여 지도를 만들어봅시다.
plot() 함수는 지도를 커스터마이징하는데 필요로 하는 데 사용할 수 있는 몇가지 파라미터(선택사항)를 입력으로 받습니다.
가장 중요한 것은 ax에 값을 설정하면, 모든 정보가 동일한 지도에 그려진다는 것입니다.
# Base Map을 지역 경계로 설정
ax = counties.plot(figsize = (10, 10), color = 'none', edgecolor = 'gainsboro', zorder = 3)
# 토지와 캠프장 위치, 도보 경로를 Base Map에 추가
wild_lands.plot(color = 'lightgreen', ax = ax)
campsites.plot(color = 'maroon', markersize = 2, ax = ax)
trails.plot(color = 'black', markersize = 1, ax = ax)<Axes: >
북동부 지역 이 캠핑 여행지에 좋은 선택지가 될 것 같습니다!
처음엔 복잡하게 느껴지겠지만, 이미 중요한 분석을 수행할 수 있을만큼 충분히 배웠을 것입니다.
Exercise1 : Your First Map를 통해 비영리 단체가 사업을 확장할 수 있는 필리핀의 외딴 지역을 직접 찾아보세요.
이 강좌에서 만드는 지도는 지구 표면을 2차원으로 묘사합니다. 하지만 아시다시피 지구는 실제로 3차원입니다. 따라서, 지구를 평면으로 렌더링하려면 Map Projection(지도 투영) 이라는 방법을 사용해야합니다.
Map Projection은 100% 정확할 수 없습니다. 각 투영법은 지구 표면을 중요한 속성은 유지하지만, 어떤 식으로든 왜곡합니다.
예를 들어,
투영된 점이 지구상의 실제 위치와 어떻게 일치하는지 보여주기 위해 좌표 참조 체계 (CRS) 를 사용합니다. 이 튜토리얼에서는 좌표 참조계에 대해 자세히 알아보고 GeoPandas에서 좌표 참조계를 사용하는 방법을 알아봅시다.
shapefile에서 GeoDataFrame을 만들 때, CRS는 이미 로드되어 있습니다.
# 가나의 지역이 포함된 GeoDataFrame 불러오기
regions = gpd.read_file('./geospatial-learn-course-data/ghana/ghana/Regions/Map_of_Regions_in_Ghana.shp')
print(regions.crs)epsg:32630이를 어떻게 해석해야하나요?
좌표 참조 시스템은 유럽 석유 측량 그룹 (EPSG) 코드에 의해 참조됩니다.
이 GeoDataFrame은 일반적으로 ‘Mercator (메르카토르)’ 투영법이라고 더 많이 불리는 EPSG 32630을 사용합니다. 이 투영법은 각도를 보존하고 (해상 항해에 유용) 면적을 약간 왜곡합니다.
그러나, CSV 파일에서 GeoDataFrame을 만들 때, CRS를 설정해야합니다. EPSG 4326은 위도 및 경도 좌표에 해당합니다.
# 가나의 의료 시설로 DataFrame을 만들기
facilities_df = pd.read_csv('./geospatial-learn-course-data/ghana/ghana/health_facilities.csv')
# DataFrame을 GeoDataFrame으로 변환하기
facilities = gpd.GeoDataFrame(facilities_df, geometry = gpd.points_from_xy(facilities_df['Longitude'], facilities_df['Latitude']))
# 좌표 참조 시스템 (CRS) EPSG 4326으로 설정
facilities.crs = {'init': 'epsg:4326'}
# GeoDataFrame의 첫 5개 행 출력
facilities.head()| Region | District | FacilityName | Type | Town | Ownership | Latitude | Longitude | geometry | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | Ashanti | Offinso North | A.M.E Zion Clinic | Clinic | Afrancho | CHAG | 7.40801 | -1.96317 | POINT (-1.96317 7.40801) |
| 1 | Ashanti | Bekwai Municipal | Abenkyiman Clinic | Clinic | Anwiankwanta | Private | 6.46312 | -1.58592 | POINT (-1.58592 6.46312) |
| 2 | Ashanti | Adansi North | Aboabo Health Centre | Health Centre | Aboabo No 2 | Government | 6.22393 | -1.34982 | POINT (-1.34982 6.22393) |
| 3 | Ashanti | Afigya-Kwabre | Aboabogya Health Centre | Health Centre | Aboabogya | Government | 6.84177 | -1.61098 | POINT (-1.61098 6.84177) |
| 4 | Ashanti | Kwabre | Aboaso Health Centre | Health Centre | Aboaso | Government | 6.84177 | -1.61098 | POINT (-1.61098 6.84177) |
위의 셀 CSV 파일에서 GeoDataFrame을 만들려면 Pandas와 GeoPandas를 모두 사용해야했습니다.
gpd.GeoDataFrame()을 사용합니다.gpd.points_from_xy() 함수는 위도 및 경도 열에서 Point 개체를 생성합니다.재투영은 CRS를 변경하는 과정을 말합니다. 이 작업은 GeoPandas에서 to_crs() 함수를 사용하여 수행됩니다.
여러 개의 GeoDataFrame을 plot할 때 모두 동일한 CRS를 사용하는 것이 중요합니다. 아래 셀에서는 facilities GeoDataFrame의 CRS를 regions의 CRS와 일치하도록 변경한 후 plot합니다.
# 지도 만들기
ax = regions.plot(figsize = (8, 8), color = 'whitesmoke', linestyle = ':', edgecolor = 'black')
facilities.to_crs(epsg = 32630).plot(markersize = 1, ax = ax)<Axes: >
to_crs() 함수는 geometry 열만 수정하고 다른 모든 열은 그대로 유지합니다.
# `Latitude`와 `Longitude` 열은 변경 X
facilities.to_crs(epsg = 32630).head()| Region | District | FacilityName | Type | Town | Ownership | Latitude | Longitude | geometry | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | Ashanti | Offinso North | A.M.E Zion Clinic | Clinic | Afrancho | CHAG | 7.40801 | -1.96317 | POINT (614422.662 818986.851) |
| 1 | Ashanti | Bekwai Municipal | Abenkyiman Clinic | Clinic | Anwiankwanta | Private | 6.46312 | -1.58592 | POINT (656373.863 714616.547) |
| 2 | Ashanti | Adansi North | Aboabo Health Centre | Health Centre | Aboabo No 2 | Government | 6.22393 | -1.34982 | POINT (682573.395 688243.477) |
| 3 | Ashanti | Afigya-Kwabre | Aboabogya Health Centre | Health Centre | Aboabogya | Government | 6.84177 | -1.61098 | POINT (653484.490 756478.812) |
| 4 | Ashanti | Kwabre | Aboaso Health Centre | Health Centre | Aboaso | Government | 6.84177 | -1.61098 | POINT (653484.490 756478.812) |
GeoPandas에서 EPSG 코드를 사용할 수 없는 경우, CRS의 ’proj4 string’을 사용하여 CRS를 변경할 수 있습니다.
예를 들어, 위도/경도 좌표로 변환하는 proj4 string은 다음과 같습니다.
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs# CRS를 EPSG 4326로 변경
regions.to_crs('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs').head()| Region | geometry | |
|---|---|---|
| 0 | Ashanti | POLYGON ((-1.30985 7.62302, -1.30786 7.62198, ... |
| 1 | Brong Ahafo | POLYGON ((-2.54567 8.76089, -2.54473 8.76071, ... |
| 2 | Central | POLYGON ((-2.06723 6.29473, -2.06658 6.29420, ... |
| 3 | Eastern | POLYGON ((-0.21751 7.21009, -0.21747 7.20993, ... |
| 4 | Greater Accra | POLYGON ((0.23456 6.10986, 0.23484 6.10974, 0.... |
첫번째 튜토리얼에서 배운 것처럼, 임의의 GeoDataFrame의 경우 geometry 열의 유형은 표시하려는 대상에 따라 달라집니다.
세 가지 유형의 기하학적 객체 모두 데이터 집합을 빠르게 분석하는데 사용할 수 있는 기본 속성이 있습니다.
예를 들어, x와 y 속성에서 각 Point의 x좌표와 y좌표를 얻을 수 있습니다.
# 각 Point에서 x좌표
facilities['geometry'].head().x0 -1.96317
1 -1.58592
2 -1.34982
3 -1.61098
4 -1.61098
dtype: float64그리고, length 속성에서 LineString의 길이를 얻을 수 있습니다.
또는, area 속성에서 Polygon의 면적을 구할 수 있습니다.
# GeoDataFrame에서 각 Polygon의 면적(km²) 계산
regions.loc[:, 'AREA'] = regions['geometry'].area / 10**6
print('Area of Ghana : {} km²'.format(regions['AREA'].sum()))
print('CRS : ', regions.crs)
regions.head()Area of Ghana : 239584.5760055668 km²
CRS : epsg:32630| Region | geometry | AREA | |
|---|---|---|---|
| 0 | Ashanti | POLYGON ((686446.075 842986.894, 686666.193 84... | 24379.017777 |
| 1 | Brong Ahafo | POLYGON ((549970.457 968447.094, 550073.003 96... | 40098.168231 |
| 2 | Central | POLYGON ((603176.584 695877.238, 603248.424 69... | 9665.626760 |
| 3 | Eastern | POLYGON ((807307.254 797910.553, 807311.908 79... | 18987.625847 |
| 4 | Greater Accra | POLYGON ((858081.638 676424.913, 858113.115 67... | 3706.511145 |
위의 셀에서, regions GeoDataFrame의 CRS가 EPSG 32630(메르카토르 투영법)으로 설정되어 있기 때문에 면적 계산은 ’Africa Albers Equal Area Conic’과 같은 등면적 투영법을 사용한 경우보다 약간 덜 정확합니다.
하지만, 이렇게 하면 가나의 면적은 약 239,584km²로 계산되며, 이는 정답과 크게 다르지 않습니다.
Exercise2 : Coordinate Reference Systems를 배운 내용을 사용하여 남미로 이주하는 조류를 추적해보세요.
이 튜토리얼에서는 folium 패키지로 대화형 지도 (Interactive Maps) 를 만드는 방법을 배웁니다. 이 과정에서, 새로운 기술을 적용하여 보스턴 범죄 데이터를 시각화합니다.
# MapBox API를 활용하여 지도 만들기
map = folium.Map(location=[36.3504119, 127.3845475], zoom_start=14,
tiles='https://api.mapbox.com/styles/v1/jw0112/cl1yi61ue002z14ojne3qrry8/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiancwMTEyIiwiYSI6ImNsaDhzcDhkMTAwcGwzY2xnenhscjZxaTMifQ.kyL3owKwnt5Ewy58TbmlGQ',
attr='Mapbox Attribution')
mapfolium.Map() 함수로 비교적 간단한 지도를 만들어보겠습니다.
# 지도 생성
m_1 = folium.Map(
location = [42.32, -71.0589],
tiles = 'openstreetmap',
zoom_start = 10)
# 지도 표현
m_1
여러 인수들을 활용하여 지도의 모양을 커스터마이징할 수 있습니다.
location은 지도의 초기 중심을 설정합니다. 여기서는 보스턴 시의 위도(42.32° N)와 경도(-71.0589° E)를 사용합니다.tiles은 지도의 스타일을 변경합니다. 여기서는 OpenStreetMap스타일을 사용합니다. 다른 옵션이 궁금하다면 이 링크에서 나열된 다른 옵션을 찾을 수 있습니다.zoom_start은 지도의 초기 줌 레벨을 설정하며, 값이 클 수록 지도가 더 가까이 확대됩니다.이제, 지도에 범죄 데이터를 추가하겠습니다!
여기서는 데이터 불러오는 것에 집중하지 않고 이미 crimes라는 Pandas DataFrame을 이미 가지고 있다고 생각하고, DataFrame의 처음 5개의 행을 보겠습니다.
# 데이터 불러오기
crimes = pd.read_csv('./geospatial-learn-course-data/crimes-in-boston/crimes-in-boston/crime.csv', encoding = 'latin-1')
# 결측치가 존재하는 행 제거
crimes.dropna(subset = ['Lat', 'Long', 'DISTRICT'], inplace = True)
# 2018년의 주요 범죄만 추출
crimes = crimes[crimes['OFFENSE_CODE_GROUP'].isin([
'Larceny', 'Auto Theft', 'Robbery', 'Larceny From Motor Vehicle',
'Residential Burglary', 'Simple Assault', 'Harassment', 'Ballistics',
'Aggravated Assault', 'Other Burglary', 'Arson', 'Commercial Burglary',
'HOME INVASION', 'Homicide', 'Criminal Harassment', 'Manslaughter'])]
crimes = crimes[crimes['YEAR'] >= 2018]
# 처음 5개 행 출력
crimes.head()| INCIDENT_NUMBER | OFFENSE_CODE | OFFENSE_CODE_GROUP | OFFENSE_DESCRIPTION | DISTRICT | REPORTING_AREA | SHOOTING | OCCURRED_ON_DATE | YEAR | MONTH | DAY_OF_WEEK | HOUR | UCR_PART | STREET | Lat | Long | Location | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | I182070945 | 619 | Larceny | LARCENY ALL OTHERS | D14 | 808 | NaN | 2018-09-02 13:00:00 | 2018 | 9 | Sunday | 13 | Part One | LINCOLN ST | 42.357791 | -71.139371 | (42.35779134, -71.13937053) |
| 6 | I182070933 | 724 | Auto Theft | AUTO THEFT | B2 | 330 | NaN | 2018-09-03 21:25:00 | 2018 | 9 | Monday | 21 | Part One | NORMANDY ST | 42.306072 | -71.082733 | (42.30607218, -71.08273260) |
| 8 | I182070931 | 301 | Robbery | ROBBERY - STREET | C6 | 177 | NaN | 2018-09-03 20:48:00 | 2018 | 9 | Monday | 20 | Part One | MASSACHUSETTS AVE | 42.331521 | -71.070853 | (42.33152148, -71.07085307) |
| 19 | I182070915 | 614 | Larceny From Motor Vehicle | LARCENY THEFT FROM MV - NON-ACCESSORY | B2 | 181 | NaN | 2018-09-02 18:00:00 | 2018 | 9 | Sunday | 18 | Part One | SHIRLEY ST | 42.325695 | -71.068168 | (42.32569490, -71.06816778) |
| 24 | I182070908 | 522 | Residential Burglary | BURGLARY - RESIDENTIAL - NO FORCE | B2 | 911 | NaN | 2018-09-03 18:38:00 | 2018 | 9 | Monday | 18 | Part One | ANNUNCIATION RD | 42.335062 | -71.093168 | (42.33506218, -71.09316781) |
지도에 넣어야 하는 데이터의 양을 줄이기 위해 (일시적으로) 낮 시간대 강도 사건으로만 한정하겠습니다.
daytime_robberies = crimes[((crimes['OFFENSE_CODE_GROUP'] == 'Robbery') & (crimes['HOUR'].isin(range(9, 18))))]
daytime_robberies.head(2)| INCIDENT_NUMBER | OFFENSE_CODE | OFFENSE_CODE_GROUP | OFFENSE_DESCRIPTION | DISTRICT | REPORTING_AREA | SHOOTING | OCCURRED_ON_DATE | YEAR | MONTH | DAY_OF_WEEK | HOUR | UCR_PART | STREET | Lat | Long | Location | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 299 | I182070598 | 311 | Robbery | ROBBERY - COMMERCIAL | A15 | 60 | NaN | 2018-09-02 17:15:00 | 2018 | 9 | Sunday | 17 | Part One | RUTHERFORD AVE | 42.371673 | -71.063264 | (42.37167264, -71.06326413) |
| 527 | I182070342 | 381 | Robbery | ROBBERY - CAR JACKING | E18 | 490 | NaN | 2018-09-01 17:05:00 | 2018 | 9 | Saturday | 17 | Part One | CUMMINS HWY | 42.276453 | -71.112980 | (42.27645319, -71.11297971) |
folium.Marker()folium.Marker() 함수로 지도에 Marker를 추가합니다. 아래의 각 Marker는 서로 다른 강도에 해당합니다.
# 지도 생성
m_2 = folium.Map(
location = [42.32, -71.0589],
tiles = 'cartodbpositron',
zoom_start = 13)
# 지도 + Marker
for idx, row in daytime_robberies.iterrows() :
Marker([row['Lat'], row['Long']]).add_to(m_2)
# 지도 표현
m_2folium.plugins.MarkerCluster()folium.plugins.MarkerCluster()를 사용하면 맵을 깔끔하게 정리할 수 있습니다. 각 Marker는 MarkerCluster 객체에 추가됩니다.
# 지도 생성
m_3 = folium.Map(
location = [42.32, -71.0589],
tiles = 'cartodbpositron',
zoom_start = 13)
# 지도 + MarkerCluster
mc = MarkerCluster()
for idx, row in daytime_robberies.iterrows() :
if not math.isnan(row['Long']) and not math.isnan(row['Lat']):
mc.add_child(Marker([row['Lat'], row['Long']]))
m_3.add_child(mc)
# Display the map
m_3버블 맵는 Marker 대신 원을 사용합니다. 각 원의 크기와 색상을 변경하여 위치와 다른 두 변수 사이의 관계를 표시할 수도 있습니다.
folium.Circle()을 사용하여 원을 반복적으로 추가하여 버블 지도를 만듭니다.
아래 셀에서 9~12시에 발생한 강도는 녹색으로 표시되고, 13~17시에 발생한 강도는 빨간색으로 표시됩니다.
# 지도 생성
m_4 = folium.Map(location = [42.32, -71.0589], tiles = 'cartodbpositron', zoom_start = 13)
def color_producer(val) :
if val <= 12 :
return 'forestgreen'
else :
return 'darkred'
# 지도 + 버블맵
for i in range(0, len(daytime_robberies)) :
Circle(
location = [daytime_robberies.iloc[i]['Lat'],
daytime_robberies.iloc[i]['Long']],
radius = 20,
color = color_producer(daytime_robberies.iloc[i]['HOUR'])).add_to(m_4)
# 지도 표현
m_4folium.Circle()은 여러 인수를 받습니다.
location은 원의 중심을 위도와 경도로 포함하는 목록입니다.radius는 원의 반지름을 설정합니다.
color_producer() 함수와 유사한 함수를 정의하여 이를 구현할 수 있습니다.color는 각 원의 색상을 설정합니다.
color_producer() 함수는 강도의 위치에 대한 시간의 효과를 시각화하는데 사용됩니다.히트맵을 만드려면 folium.plugins.HeatMap()을 사용합니다. 이는 도시 내 여러 지역의 범죄 밀도를 보여주며, 빨간색 영역은 상대적으로 범죄 발생이 더 많습니다.
대도시에서 예상할 수 있듯이, 대부분의 범죄는 도심 근처에서 발생합니다.
# 지도 생성
m_5 = folium.Map(location = [42.32, -71.0589], tiles = 'cartodbpositron', zoom_start = 12)
# 지도 + 히트맵
HeatMap(data = crimes[['Lat', 'Long']], radius = 10).add_to(m_5)
# 지도 표현
m_5folium.plugins.HeatMap()은 여러 인수를 받습니다.
data는 시각화하고자 하는 위치를 포함하는 DataFrame입니다.radius는 히트맵의 부드러움 정도를 조정합니다. 값이 클수록 히트맵의 부드러움 정도가 커집니다.(즉, 간격이 줄어듭니다.)경찰 관할 구역별로 범죄가 어떻게 다른지 이해하기 위해 단계구분도를 만들어 보겠습니다.
첫번째 단계로, 각 구역에서 서로 다른 행이 할당되고 geometry 열에 지리적 경계가 포함된 GeoDataFrame을 만듭니다.
# 보스턴 경찰 관할 구역 지리적 경계가 포함된 GeoDataFrame
districts_full = gpd.read_file('./geospatial-learn-course-data/Police_Districts/Police_Districts/Police_Districts.shp')
districts = districts_full[['DISTRICT', 'geometry']].set_index('DISTRICT')
districts.head()| geometry | |
|---|---|
| DISTRICT | |
| A15 | MULTIPOLYGON (((-71.07416 42.39051, -71.07415 ... |
| A7 | MULTIPOLYGON (((-70.99644 42.39557, -70.99644 ... |
| A1 | POLYGON ((-71.05200 42.36884, -71.05169 42.368... |
| C6 | POLYGON ((-71.04406 42.35403, -71.04412 42.353... |
| D4 | POLYGON ((-71.07416 42.35724, -71.07359 42.357... |
각 구역별 범죄 발생 건수를 보여주는 plot_dict라는 Pandas Series를 만듭니다.
# 각 구역별 범죄 발생 건수
plot_dict = crimes['DISTRICT'].value_counts()
plot_dict.head()D4 2885
B2 2231
A1 2130
C11 1899
B3 1421
Name: DISTRICT, dtype: int64지리적 경계를 적절한 색상과 일치시키는 방법을 알기 때문에, plot_dict과 districts가 동일한 DISTRICT라는 인덱스를 갖는 것이 매우 중요합니다.
folium.Choropleth() 함수를 사용하여 단계구분도를 만들 수 있습니다. 아래 지도가 렌더링되지 않는 경우 다른 웹 브라우저에서 페이지를 확인해보세요.
# 지도 생성
m_6 = folium.Map(location = [42.32, -71.0589], tiles = 'cartodbpositron', zoom_start = 12)
# 지도 + 단계구분도
Choropleth(geo_data = districts.__geo_interface__,
data = plot_dict,
key_on = 'feature.id',
fill_color = 'YlGnBu',
legend_name = 'Major criminal incidents (Jan-Aug 2018)'
).add_to(m_6)
# 지도 표현
m_6folium.Choropleth()는 여러 인수를 받습니다.
geo_data는 각 지리적 영역의 경계를 포함하는 GeoJSON FeatureCollection입니다.
districts GeoDataFrame을 __geo_interface__ 속성을 사용하여 GeoJSON FeatureCollection으로 변환합니다.data는 각 지리적 영역의 색상을 지정하는데 사용할 값이 포함된 Pandas Series입니다.key_on은 항상 feature.id로 설정됩니다.
geo_data에 사용되는 GeoDataFrame과 data에 사용되는 Pandas Series가 동일한 인덱스를 가지고 있다는 사실을 나타냅니다. 자세한 내용을 이해하려면 GeoJSON Feature Collection의 구조를 좀 더 자세히 살펴봐야합니다. (‘feature’ 키에 해당하는 값은 list이고, 각 항목은 ‘id’ 키가 포함된 dictonary입니다.)fill_color는 색상 스케일을 설정합니다.legend_name은 맵 오른쪽 상단에 있는 범례에 레이블을 지정합니다.Exercise3 : Interactive Maps을 통해 일본에서 추가적으로 내진 보강이 필요한 지역을 확인해보세요.
이 튜토리얼에서는 지리공간 데이터에 대한 (지오코딩 (Geocoding) 과 테이블 조인 (Table Joins)) 두가지 조작에 대해 알아봅시다.
지오코딩은 장소 이름이나 주소를 지도 상의 위치로 변환하는 과정입니다. 예를 들어, Google Maps, Bing Maps, Baidu Maps에서 랜드마크 설명을 기반으로 지리적 위치를 찾아본 적이 있다면 지오코딩을 사용해 본 적이 있을 것입니다!
모든 지오코딩에는 geopy가 사용됩니다.
from geopy.geocoders import Nominatim위의 셀에서 Nominatim은 위치를 생성하는데 사용되는 지오코딩 소프트웨어를 나타냅니다.
지오코더를 인스턴트화하는 것으로 시작합니다. 그런 다음, 이름이나 주소를 Python 문자열로 적용하기만 하면 됩니다. (이 경우, 기자의 대피라미드라고 알려진 "Pyramid of Khufu"를 제공합니다.)
지오코딩이 성공하면 두 가지 중요한 속성을 가진 geopy.location.Location 객체를 반환합니다.
point 속성에는 (위도, 경도) 위치가 포함합니다.address 속성에는 전체 주소를 포함합니다.geolocator = Nominatim(user_agent = 'kaggle_learn')
location = geolocator.geocode('Pyramid of Khufu')
print(location.point)
print(location.address)29 58m 44.976s N, 31 8m 3.17625s E
هرم خوفو, شارع ابو الهول السياحي, نزلة البطران, الجيزة, 12125, مصرpoint 속성의 값은 geopy.point.Point 객체이며 각각 latitude와 longitude 속성에서 위도와 경도를 가져올 수 있습니다.
point = location.point
print('Latitude : ', point.latitude)
print('Longitude : ', point.longitude)Latitude : 29.97916
Longitude : 31.134215625236113종종 다양한 주소를 지오코딩해야 하는 경우도 있습니다. 예를 들어, 유럽에 있는 상위 100개 대학의 위치를 얻고 싶다고 가정해 보겠습니다.
universities = pd.read_csv('./geospatial-learn-course-data/top_universities.csv')
universities.head()| Name | |
|---|---|
| 0 | University of Oxford |
| 1 | University of Cambridge |
| 2 | Imperial College London |
| 3 | ETH Zurich |
| 4 | UCL |
그런 다음 lambda 함수를 사용해서 DataFrame의 모든 행에 지오코딩 적용할 수 있습니다. (지오코딩이 실패할 경우를 대비하여 try / except 문을 사용합니다.)
def my_geocoder(row) :
try :
point = geolocator.geocode(row).point
return pd.Series({'Latitude' : point.latitude, 'Longitude' : point.longitude})
except :
return None
universities[['Latitude', 'Longitude']] = universities.apply(lambda x : my_geocoder(x['Name']), axis = 1)
print('{}% of addresses were geocoded!'.format(
(1 - sum(np.isnan(universities['Latitude'])) / len(universities)) * 100))
# 지오코딩 실패한 행 제거
universities = universities.loc[~np.isnan(universities['Latitude'])]
universities = gpd.GeoDataFrame(
universities, geometry = gpd.points_from_xy(universities.Longitude, universities.Latitude))
universities.crs = {'init' : 'epsg:4326'}
universities.head()91.0% of addresses were geocoded!| Name | Latitude | Longitude | geometry | |
|---|---|---|---|---|
| 0 | University of Oxford | 51.759037 | -1.252430 | POINT (-1.25243 51.75904) |
| 1 | University of Cambridge | 52.200623 | 0.110474 | POINT (0.11047 52.20062) |
| 2 | Imperial College London | 51.498959 | -0.175641 | POINT (-0.17564 51.49896) |
| 3 | ETH Zurich | 47.562772 | 7.580947 | POINT (7.58095 47.56277) |
| 4 | UCL | 51.521785 | -0.135151 | POINT (-0.13515 51.52179) |
다음으로, 지오코더가 반환한 모든 위치를 시각화합니다. 몇몇 위치는 유럽에 있지 않기 때문에 확실히 부정확하다는 것을 알 수 있습니다!
# 지도 생성
m = folium.Map(location = [54, 15], tiles = 'openstreetmap', zoom_start = 2)
# 지도 + Marker
for idx, row in universities.iterrows() :
Marker([row['Latitude'], row['Longitude']], popup = row['Name']).add_to(m)
# 지도 표현
m이제 주제를 전환하여 서로 다른 원본의 데이터를 결합하는 방법에 대해 알아보겠습니다.
이미 pd.DataFrame.join()을 사용하여 여러 DataFrame의 정보를 동일한 인덱스로 결합하는 방법을 알고 있습니다. 인덱스에서 일치하는 값을 단순화하여 데이터를 조인하는 방법을 특성 조인 (Attribute Join)이라고 합니다.
GeoDataFrame의 특성 조인을 수행할 때는 gpd.GeoDataFrame.merge()을 사용하는 것이 가장 좋습니다. 이를 설명하기 위해, 유럽의 모든 국가에 대한 경계가 포함된 europe_boundaries GeoDataFrame으로 작업해보겠습니다.
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
europe = world.loc[world.continent == 'Europe'].reset_index(drop = True)
europe_stats = europe[['name', 'pop_est', 'gdp_md_est']]
europe_boundaries = europe[['name', 'geometry']]
europe_boundaries.head()| name | geometry | |
|---|---|---|
| 0 | Russia | MULTIPOLYGON (((180.00000 71.51571, 180.00000 ... |
| 1 | Norway | MULTIPOLYGON (((15.14282 79.67431, 15.52255 80... |
| 2 | France | MULTIPOLYGON (((-51.65780 4.15623, -52.24934 3... |
| 3 | Sweden | POLYGON ((11.02737 58.85615, 11.46827 59.43239... |
| 4 | Belarus | POLYGON ((28.17671 56.16913, 29.22951 55.91834... |
각 국가의 예상 인구와 국내총생산(GDP)을 포함하는 DataFrame europe_stats와 결합하겠습니다.
europe_stats.head()| name | pop_est | gdp_md_est | |
|---|---|---|---|
| 0 | Russia | 144373535.0 | 1699876 |
| 1 | Norway | 5347896.0 | 403336 |
| 2 | France | 67059887.0 | 2715518 |
| 3 | Sweden | 10285453.0 | 530883 |
| 4 | Belarus | 9466856.0 | 63080 |
아래 셀에서 특성 조인을 수행합니다. on 인수는 europe_boundaries의 행을 europe_stats의 행과 일치시키는데 사용되는 열 이름으로 설정됩니다.
# 특성 조인을 사용하여 유럽 국가에 대한 데이터 병합
europe = europe_boundaries.merge(europe_stats, on = 'name')
europe.head()| name | geometry | pop_est | gdp_md_est | |
|---|---|---|---|---|
| 0 | Russia | MULTIPOLYGON (((180.00000 71.51571, 180.00000 ... | 144373535.0 | 1699876 |
| 1 | Norway | MULTIPOLYGON (((15.14282 79.67431, 15.52255 80... | 5347896.0 | 403336 |
| 2 | France | MULTIPOLYGON (((-51.65780 4.15623, -52.24934 3... | 67059887.0 | 2715518 |
| 3 | Sweden | POLYGON ((11.02737 58.85615, 11.46827 59.43239... | 10285453.0 | 530883 |
| 4 | Belarus | POLYGON ((28.17671 56.16913, 29.22951 55.91834... | 9466856.0 | 63080 |
또 다른 조인 유형은 공간 조인(Spatial Join)입니다. 공간 조인을 사용하면 geometry 열에 있는 개체 간의 공간 관계를 기반으로 GeoDataFrame을 결합합니다. 예를 들어, 유럽 대학의 지오코딩된 주소가 포함된 GeoDataFrame universities가 이미 있습니다.
그런 다음 공간 조인을 사용하여 각 대학을 해당 국가에 일치시킬 수 있습니다. 이 작업은 gpd.sjoin()을 사용하여 수행합니다.
# 공간 조인을 사용하여 대학을 유럽의 국가와 일치
european_universities = gpd.sjoin(universities, europe)
# 결과 검토
print('We located {} universities.'.format(len(universities)))
print('Only {} of the universities were located in Europe (in {} different countries).'.format(
len(european_universities), len(european_universities.name.unique())))
european_universities.head()We located 91 universities.
Only 86 of the universities were located in Europe (in 14 different countries).| Name | Latitude | Longitude | geometry | index_right | name | pop_est | gdp_md_est | |
|---|---|---|---|---|---|---|---|---|
| 0 | University of Oxford | 51.759037 | -1.252430 | POINT (-1.25243 51.75904) | 28 | United Kingdom | 66834405.0 | 2829108 |
| 1 | University of Cambridge | 52.200623 | 0.110474 | POINT (0.11047 52.20062) | 28 | United Kingdom | 66834405.0 | 2829108 |
| 2 | Imperial College London | 51.498959 | -0.175641 | POINT (-0.17564 51.49896) | 28 | United Kingdom | 66834405.0 | 2829108 |
| 4 | UCL | 51.521785 | -0.135151 | POINT (-0.13515 51.52179) | 28 | United Kingdom | 66834405.0 | 2829108 |
| 5 | London School of Economics and Political Science | 51.514211 | -0.116808 | POINT (-0.11681 51.51421) | 28 | United Kingdom | 66834405.0 | 2829108 |
위의 공간 조인은 두 GeoDataFrame의 geometry 열을 살펴봅니다. universities GeoDataFrame의 Point 객체가 europe DataFrame의 Polygon 객체와 교차하는 경우, 해당 행이 결합되어 european_universities DataFrame의 단일 행으로 추가됩니다. 그렇지 않으면, 일치하는 대학이 없는 국가(및 일치하는 국가가 없는 대학)는 결과에서 제거됩니다.
gpd.sjoin() 함수는 how 및 op 인수를 통해 다양한 조인 유형에 맞게 지정할 수 있습니다. 예를 들어, how = left(또는 how = right)를 설정하여 SQL Left (또는 Right) 조인과 동등한 작업을 수행할 수 있습니다. 이 강좌에서는 자세히 설명하지 않겠지만, 이 문서에서 자세한 내용을 확인할 수 있습니다.
Exercise 4 : Manipulating Geospatial Data 에서 지오코딩 및 테이블 조인 사용을 사용하여 스타벅스 리저브 로스터리에 적합한 위치를 파악해봅시다.
이 튜토리얼에서는 근접성 분석 을 위한 몇가지 기법을 살펴봅니다. 특히 다음과 같은 작업을 수행하는 방법을 배우게 됩니다.
미국 펜실베니아주 필라델피아의 독성 화학물질 배출을 추적하는 미국 환경보호청 (EPA)의 데이터를 사용하여 작업하겠습니다.
releases = gpd.read_file('./geospatial-learn-course-data/toxic_release_pennsylvania/toxic_release_pennsylvania/toxic_release_pennsylvania.shp')
releases.head()| YEAR | CITY | COUNTY | ST | LATITUDE | LONGITUDE | CHEMICAL | UNIT_OF_ME | TOTAL_RELE | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2016 | PHILADELPHIA | PHILADELPHIA | PA | 40.005901 | -75.072103 | FORMIC ACID | Pounds | 0.160 | POINT (2718560.227 256380.179) |
| 1 | 2016 | PHILADELPHIA | PHILADELPHIA | PA | 39.920120 | -75.146410 | ETHYLENE GLYCOL | Pounds | 13353.480 | POINT (2698674.606 224522.905) |
| 2 | 2016 | PHILADELPHIA | PHILADELPHIA | PA | 40.023880 | -75.220450 | CERTAIN GLYCOL ETHERS | Pounds | 104.135 | POINT (2676833.394 261701.856) |
| 3 | 2016 | PHILADELPHIA | PHILADELPHIA | PA | 39.913540 | -75.198890 | LEAD COMPOUNDS | Pounds | 1730.280 | POINT (2684030.004 221697.388) |
| 4 | 2016 | PHILADELPHIA | PHILADELPHIA | PA | 39.913540 | -75.198890 | BENZENE | Pounds | 39863.290 | POINT (2684030.004 221697.388) |
또한, 같은 도시의 대기질 모니터링 스테이션에서 측정한 수치가 포함된 데이터로 작업하겠습니다.
stations = gpd.read_file('./geospatial-learn-course-data/PhillyHealth_Air_Monitoring_Stations/PhillyHealth_Air_Monitoring_Stations/PhillyHealth_Air_Monitoring_Stations.shp')
stations.head()| SITE_NAME | ADDRESS | BLACK_CARB | ULTRAFINE_ | CO | SO2 | OZONE | NO2 | NOY_NO | PM10 | ... | PAMS_VOC | TSP_11101 | TSP_METALS | TSP_LEAD | TOXICS_TO1 | MET | COMMUNITY_ | LATITUDE | LONGITUDE | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | LAB | 1501 East Lycoming Avenue | N | N | Y | N | Y | Y | Y | N | ... | Y | N | Y | N | y | N | N | 40.008606 | -75.097624 | POINT (2711384.641 257149.310) |
| 1 | ROX | Eva and Dearnley Streets | N | N | N | N | N | N | N | N | ... | N | N | Y | N | Y | N | N | 40.050461 | -75.236966 | POINT (2671934.290 271248.900) |
| 2 | NEA | Grant Avenue and Ashton Street | N | N | N | N | Y | N | N | N | ... | N | N | N | N | N | Y | N | 40.072073 | -75.013128 | POINT (2734326.638 280980.247) |
| 3 | CHS | 500 South Broad Street | N | N | N | N | N | N | N | N | ... | N | N | Y | N | Y | N | N | 39.944510 | -75.165442 | POINT (2693078.580 233247.101) |
| 4 | NEW | 2861 Lewis Street | N | N | Y | Y | Y | N | Y | Y | ... | N | Y | N | Y | N | Y | N | 39.991688 | -75.080378 | POINT (2716399.773 251134.976) |
5 rows × 24 columns
서로 다른 두 GeoDataFrame의 Point 간의 거리를 측정하려면, 먼저 두 GeoDataFrame이 동일한 좌표 참조 시스템(CRS)을 사용하는지 확인해야 합니다.
다행히도 여기서는 둘 다 EPSG:2272를 사용합니다.
print(stations.crs)
print(releases.crs)epsg:2272
epsg:2272또한, CRS가 어떤 단위(meter, feet 또는 다른 단위)를 사용하는지 확인합니다. 이 경우 EPSG:2272는 feet 단위를 사용합니다. (원한다면 여기에서 확인할 수 있습니다.)
GeoPandas에서 거리를 계산하는 것은 비교적 간단합니다. 아래 셀은 recent_release의 비교적 최근 릴리즈 사건과 stations GeoDataFrame의 모든 스테이션 사이의 거리(피트)를 계산합니다.
# 특정 릴리즈 사건 1개 선택
recent_release = releases.iloc[360]
# 릴리즈에서 각 스테이션까지의 거리 측정
distances = stations['geometry'].distance(recent_release['geometry'])
distances0 44778.509761
1 51006.456589
2 77744.509207
3 14672.170878
4 43753.554393
5 4711.658655
6 23197.430858
7 12072.823097
8 79081.825506
9 3780.623591
10 27577.474903
11 19818.381002
dtype: float64계산된 거리를 사용하여, 각 스테이션까지의 평균 거리와 같은 통계를 얻을 수 있습니다.
print('Mean distance to monitoring stations: {} feet'.format(distances.mean()))Mean distance to monitoring stations: 33516.28487007786 feet또한, 가장 가까운 모니터링 스테이션을 찾을 수도 있습니다.
print('Closest monitoring station ({} feet):'.format(distances.min()))
print(stations.iloc[distances.idxmin()][['ADDRESS', 'LATITUDE', 'LONGITUDE']])Closest monitoring station (3780.623590556444 feet):
ADDRESS 3100 Penrose Ferry Road
LATITUDE 39.91279
LONGITUDE -75.185448
Name: 9, dtype: object지도에서 한 점으로부터 일정 반경 떨어진 모든 점을 이해하려면 버퍼를 만드는 것이 가장 간단한 방법입니다.
아레 셀은 12개의 서로 다른 Polygon 객체를 포함하는 GeoSeries two_mile_buffer를 만듭니다. 각 다각형은 서로 다른 대기 모니터링 스테이션 주변의 2 miles(또는, 2*5280 feet)의 버퍼입니다.
two_mile_buffer = stations['geometry'].buffer(2*5280)
two_mile_buffer.head()0 POLYGON ((2721944.641 257149.310, 2721893.792 ...
1 POLYGON ((2682494.290 271248.900, 2682443.441 ...
2 POLYGON ((2744886.638 280980.247, 2744835.789 ...
3 POLYGON ((2703638.580 233247.101, 2703587.731 ...
4 POLYGON ((2726959.773 251134.976, 2726908.924 ...
dtype: geometry각 Polygon을 지도에 그리기 위해 folium.GeoJson()을 사용합니다. folium에는 위도와 경도의 좌표가 필요하므로 시각화 전에 CRS를 EPSG:4326으로 변환해야합니다.
# 최근 릴리즈 사건과 모니터링 스테이션을 포함하여 지도를 생성
m = folium.Map(location = [39.9526, -75.1652], zoom_start = 11)
HeatMap(data = releases[['LATITUDE', 'LONGITUDE']], radius = 15).add_to(m)
for idx, row in stations.iterrows() :
Marker([row['LATITUDE'], row['LONGITUDE']]).add_to(m)
# 지도 + 버퍼
GeoJson(two_mile_buffer.to_crs(epsg = 4326)).add_to(m)
# 지도 표현
m이제 모든 모니터링 스테이션에서 2마일 이내에 독성 물질 방출이 발생했는지 테스트하려면 각 Polygon에 대해 12가지 테스트를 실행하여 해당 지점이 포함되어 있는지 개별적으로 확인할 수 있습니다.
하지만 더 효율적인 방법은 먼저 모든 Polygon을 MultiPolygon 객체로 축소하는 것입니다. 이 작업은 unary_union 속성을 사용하여 수행합니다.
# Polygon 그룹을 단일 MultiPolygon 객체로 축소
my_union = two_mile_buffer.geometry.unary_union
print('Type : ', type(my_union))
# MultiPolygon 객체 시각화
my_unionType : <class 'shapely.geometry.multipolygon.MultiPolygon'>
MultiPolygon에 Point가 포함되어 있는지 확인하기 위해 contains() 함수를 사용합니다. 튜토리얼 앞부분의 릴리즈 사건을 사용하겠습니다.
가장 가까운 모니터링 스테이션까지의 거리가 약 3,781피트라는 것을 알고 있습니다.
하지만 모든 릴리즈가 대기 모니터링 스테이션에서 2마일 이내에 발생한 것은 아닙니다!
# 가장 가까운 스테이션이 2마일 미만인 경우
print(my_union.contains(releases.iloc[360].geometry))
# 가장 가까운 스테이션이 2마일 이상 떨어진 경우
print(my_union.contains(releases.iloc[358].geometry))True
FalseExercise 5 : Proximity Analysis 에서 뉴욕 시의 병원 서비스 적용 범위를 알아보세요.